# Load packages
library(here)
library(dplyr)
library(readr)
library(rstan)
library(bayesrules)
library(tidyverse)
library(bayesplot)
library(rstanarm)
library(janitor)
library(tidybayes)
library(broom.mixed)
library(here)
library(sf)
library(tidycensus)
library(openxlsx)
library(s2)
# themes
theme_set(theme_minimal())
vari_names <- read_csv(here("clean_data", "nyc_names.csv"))
nyc_clean <- st_read(here("clean_data", "nyc_data.shp"), crs = 4269) 
## Reading layer `nyc_data' from data source 
##   `/Users/freddy/Documents/GitHub/454/clean_data/nyc_data.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 242 features and 22 fields (with 1 geometry empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25609 ymin: 40.4961 xmax: -73.70002 ymax: 40.91758
## Geodetic CRS:  NAD83
colnames(nyc_clean) <- colnames(vari_names)
subway_stations <- st_read(here("ethnic","Data","stations", "geo_export_85568705-efba-4456-bdc0-3d70ff2cf8e5.shp")) %>%
  st_transform(., 4269)
## Reading layer `geo_export_85568705-efba-4456-bdc0-3d70ff2cf8e5' from data source 
##   `/Users/freddy/Documents/GitHub/454/ethnic/Data/stations/geo_export_85568705-efba-4456-bdc0-3d70ff2cf8e5.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 473 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.03088 ymin: 40.57603 xmax: -73.7554 ymax: 40.90313
## Geodetic CRS:  WGS84(DD)
bus_stations <- st_read(here("ethnic","Data","bus", "bus_stops_nyc_may2020.shp")) %>%
  st_transform(., 4269)
## Reading layer `bus_stops_nyc_may2020' from data source 
##   `/Users/freddy/Documents/GitHub/454/ethnic/Data/bus/bus_stops_nyc_may2020.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 14663 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 914169.1 ymin: 122626.8 xmax: 1066982 ymax: 271696.8
## Projected CRS: NAD83 / New York Long Island (ftUS)
transit_points <- read_csv(here("transit","ridership_points.csv"))%>%
  separate(Position, into=c("Point", "longitude", "latitude"), " ") %>%
  mutate(latitude = str_remove_all(latitude, "[)]"),
         longitude = str_remove_all(longitude, "[()]"),
         ) %>%
  dplyr::select(-c(Point)) %>% 
  mutate(latitude = as.numeric(latitude),
         longitude = as.numeric(longitude)) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4269)
#plot locations over map
subway_loc <- ggplot() +
  geom_sf(data = nyc_clean, fill = "#EBF6FF", color = "#D48DD8", size = 0.15, alpha = .8) +
  geom_sf(data = subway_stations, color="#3F123C", size=1) + 
  coord_sf(datum = st_crs(subway_stations)) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Subway Stop Locations \nin NYC")+ 
    theme(#panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 30, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

bus_loc <- ggplot() +
  geom_sf(data = nyc_clean, fill = "#EBF6FF", color = "#D48DD8", size = 0.15, alpha = .8) +
  geom_sf(data = bus_stations, color="#3F123C", size=.5, alpha=.5) + 
  coord_sf(datum = st_crs(subway_stations)) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Bus Stop Locations \nin NYC")+ 
    theme(#panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 30, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))


stops <- nyc_clean %>%
  ggplot() +
  geom_sf(aes(fill = sub_count), color = "#8f98aa") +
  scale_fill_gradient(low= "lavender", high = "maroon",
                      guide = guide_legend(title = "Number of Subway Stops") ,na.value="#D6D6D6") +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Subway Stop Counts \nin NYC")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 30, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

bus_stops <- nyc_clean %>%
  ggplot() +
  geom_sf(aes(fill = bus_count), color = "#8f98aa") +
  scale_fill_gradient(low= "lavender", high = "maroon",
                      guide = guide_legend(title = "Number of Bus Stops") ,na.value="#D6D6D6") +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Bus Stop Counts \nin NYC")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 30, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))


ridership <- nyc_clean%>%
  ggplot() +
  geom_sf(aes(fill = log2(mean_ridership)), color = "#8f98aa") +
  scale_fill_gradient(low= "lavender", high = "maroon",
                      guide = guide_legend(title = "Log2 Mean Ridership") ,na.value="#D6D6D6") +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Mean (Log2) Subway Turnstile \nRidership in 2018 \nfor NYC")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 30, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))


access <- nyc_clean %>%
  ggplot() +
  geom_sf(aes(fill = perc_covered_by_transit), color = "#8f98aa") +
  scale_fill_gradient(low= "lavender", high = "maroon",
                      guide = guide_legend(title = "Percent of Neighborhood within \n walking distance of a transit station"), na.value="#D6D6D6") +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Transit Accessibility \nin NYC")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 30, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

`

red <- ggplot(nyc_clean) +
  geom_sf(aes(fill = below_poverty_line_count), color = "#8f98aa") +
  scale_fill_gradient(low = "#FCF5EE", high = "#E13728", guide = guide_legend(title = "Number Below \nPoverty Line")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Impoverishement")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

yellow <- ggplot(nyc_clean) +
  geom_sf(aes(fill = mean_income), color = "#8f98aa") +
  scale_fill_gradient(low = "#FCF5EE", high = "#F3D24E", guide = guide_legend(title = "Mean Income")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Mean Income")+
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12),
          legend.text = element_text(size = 12)) +
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

teal <- ggplot(nyc_clean) +
  geom_sf(aes(fill = mean_rent), color = "#8f98aa") +
  scale_fill_gradient(low = "#FCF5EE", high = "#2DBDC7", guide = guide_legend(title = "Dollars")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Mean Rent")+
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12),
          legend.text = element_text(size = 12)) +
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

purple <- ggplot(nyc_clean) +
  geom_sf(aes(fill = eviction_count), color = "#8f98aa")+
  scale_fill_gradient(low = "#FCF5EE", high = "#7826C0", guide = guide_legend(title = "Number of Evictions")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Evictions")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

orange <- ggplot(nyc_clean) +
  geom_sf(aes(fill = unemployment_count), color = "#8f98aa")+
  scale_fill_gradient(low = "#FCF5EE", high = "#FC9228", guide = guide_legend(title = "Number on \nUnemployment")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Unemployment")+
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12),
          legend.text = element_text(size = 12)) +
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

green <- ggplot(nyc_clean) +
  geom_sf(aes(fill = store_count), color = "#8f98aa")+
  scale_fill_gradient(low = "#FCF5EE", high = "#326902", guide = guide_legend(title = "Number of Stores")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Retail Food Stores")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

blue <- ggplot(nyc_clean) +
  geom_sf(aes(fill = school_count), color = "#8f98aa")+
  scale_fill_gradient(low = "#FCF5EE", high = "#5372C4", 
                      guide = guide_legend(title = "Number of Schools")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Number of Schools")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

pink <- ggplot(nyc_clean) +
  geom_sf(aes(fill = total_pop), color = "#8f98aa")+
  scale_fill_gradient(low = "#FCF5EE", high = "#F450E1", guide = guide_legend(title = "Number of People")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Population")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))


brown <- ggplot(nyc_clean) +
  geom_sf(aes(fill = uninsured_count), color = "#8f98aa")+
  scale_fill_gradient(low = "#F8E3DD", high = "#6A4D39", guide = guide_legend(title = "Number of People \n without Insurance Coverage")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Insurance Coverage")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 26, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))
white <- ggplot(nyc_clean) +
  geom_sf(aes(fill = white_count), color = "#8f98aa") +
  scale_fill_gradient(low = "#FCF5EE", high = "#7B435B", guide = guide_legend(title = "Number White")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("White Population")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 24, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

black <- ggplot(nyc_clean) +
  geom_sf(aes(fill = black_count), color = "#8f98aa") +
  scale_fill_gradient(low = "#FCF5EE", high = "#F25F5C", guide = guide_legend(title = "Number Black")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Black Population")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 24, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

asian <- ggplot(nyc_clean) +
  geom_sf(aes(fill = asian_count), color = "#8f98aa") +
  scale_fill_gradient(low = "#FCF5EE", high = "#717EC3", guide = guide_legend(title = "Number Asian")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Asian Population")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 24, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

latinx <- ggplot(nyc_clean) +
  geom_sf(aes(fill = latinx_count), color = "#8f98aa")+
  scale_fill_gradient(low = "#FCF5EE", high = "#FC9A38", guide = guide_legend(title = "Number Latinx")) +
  theme_minimal() +
  theme(panel.grid.major = element_line("transparent"),
        axis.text = element_blank()) +
  ggtitle("Latinx Population")+ 
    theme(panel.grid.major = element_line("transparent"),
          plot.title = element_text(size = 24, face = "bold"),
          legend.title = element_text(size = 12), 
          legend.text = element_text(size = 12)) + 
    guides(shape = guide_legend(override.aes = list(size = 8)),
           color = guide_legend(override.aes = list(size = 8)))

Data Introduction

All the data used in this project are from two major sources: the Tidycensus package and NYC Open Data.

Tidycensus is an R package interface, developed by Kyle Walker and Matt Herman, that enables easy access to the US Census Bureau’s data APIs and returns Tidyverse-ready data frames from various major US Census Bureau datasets. Our demographic and socioeconomic data are drawn from the American Community Survey results found in Tidycensus package. A summary of our ACS data variables is below:

  • total_pop: Total Population by Neighborhood

  • mean_income: Mean Income by Neighborhood

  • below_poverty_line_count: Number of People Living Below the 100% Poverty Line by Neighborhood

  • mean_rent: Mean Rent by Neighborhood

  • unemployment_count: Number of People on Unemployment by Neighborhood

  • latinx_count: Number of Latinx People by Neighborhood

  • white_count: Number of White People by Neighborhood

  • black_count: Number of Black People by Neighborhood

  • native_count: Number of Native People by Neighborhood

  • asian_count: Number of Asian People by Neighborhood

  • naturalized_citizen_count: Number of Naturalized Citizens by Neighborhood

  • noncitizen_count: Number of Foreign Born People by Neighborhood

  • uninsured_count: Number of Uninsured Citizens of any Age by Neighborhood

For remaining predictors, we used NYC Open Data’s portal to identify specific predictors. In particular, we used geotagged locations of Subway Stops, Bus Stops, Grocery Stores, Schools, and Eviction Sites from the Departments of Transportation, Health, Education, and Housing to calculate neighborhood-specific variables described below:

  • school_count: Number of Public Schools by Neighborhood
  • eviction_count: Number of Evictions by Neighborhood
  • store_count: Number of Grocery Stores and Food Vendors by Neighborhood
  • sub_count: Number of Subway Stations by Neighborhood
  • bus_count: Number of Bus Stations by Neighborhood
  • perc_covered_by_transit: Percent of Neighborhood within walking distance (.25 miles) of any transit stop.

Lastly, we acquired subway ridership from Metropolitan Transportation Authority’s turnstile data for the week of September 7, 2019. For each station, entry/exit of each turnstile is recorded. Then, we aggregated this information by taking the station-specific average of subway ridership across the 7 days in the week. Finally, we geotagged each listed station, then took the mean of ridership at all subway stations in each neighborhood.

Data Summaries

Our data has 242 observations of 23 variables. Below is a preview of our dataset with colnames attached.

library(kableExtra)
dim(nyc_clean)
## [1] 242  23
kable(head(nyc_clean)) %>% kable_styling() 
nta_id total_pop mean_income below_poverty_line_count below_poverty_line_and_50_count mean_rent unemployment_count latinx_count white_count black_count native_count asian_count naturalized_citizen_count noncitizen_count uninsured_count school_count eviction_count store_count sub_count bus_count mean_ridership perc_covered_by_transit geometry
BK0101 26308 98338.67 2397 1289 2062.667 582 3284 20526 482 40 1052 3777 3129 1797 6 68 71 2 53 9410.500 28.36395 MULTIPOLYGON (((-73.94074 4…
BK0102 57774 101238.92 9120 4474 2330.077 1710 18227 32237 1460 0 4008 6802 7746 3725 12 204 129 2 97 26603.000 55.26492 MULTIPOLYGON (((-73.96355 4…
BK0103 36891 30309.25 18285 5970 1194.875 457 3351 31799 1288 20 194 3548 1012 711 6 45 58 3 35 6348.667 87.71022 MULTIPOLYGON (((-73.96762 4…
BK0104 41861 78746.25 8406 3876 1801.667 1678 13602 17682 3960 210 5515 5770 5325 3040 15 157 117 6 77 8006.000 63.00854 MULTIPOLYGON (((-73.95083 4…
BK0201 23758 140543.00 1504 585 2275.833 579 1517 17643 1330 44 2041 2082 1448 852 1 25 30 2 18 5275.000 98.46067 MULTIPOLYGON (((-73.99066 4…
BK0202 24603 132850.00 3776 970 2413.875 1195 3772 13288 3369 0 2893 2151 2728 1244 17 111 73 8 104 19444.000 159.13596 MULTIPOLYGON (((-73.99327 4…

Below is a summary of each variable’s distribution.

summary(nyc_clean) 
##     nta_id            total_pop      mean_income     below_poverty_line_count
##  Length:242         Min.   :    0   Min.   : 23149   Min.   :    0           
##  Class :character   1st Qu.:16412   1st Qu.: 51287   1st Qu.: 1509           
##  Mode  :character   Median :30450   Median : 67988   Median : 3897           
##                     Mean   :31380   Mean   : 72767   Mean   : 5620           
##                     3rd Qu.:45786   3rd Qu.: 87604   3rd Qu.: 8383           
##                     Max.   :97786   Max.   :211822   Max.   :28755           
##                                     NA's   :44                               
##  below_poverty_line_and_50_count   mean_rent    unemployment_count
##  Min.   :    0.0                 Min.   : 792   Min.   :   0.0    
##  1st Qu.:  764.5                 1st Qu.:1308   1st Qu.: 352.8    
##  Median : 2302.0                 Median :1500   Median : 890.0    
##  Mean   : 2941.7                 Mean   :1583   Mean   :1032.1    
##  3rd Qu.: 4505.5                 3rd Qu.:1709   3rd Qu.:1490.0    
##  Max.   :13934.0                 Max.   :3279   Max.   :4697.0    
##                                  NA's   :44                       
##   latinx_count    white_count     black_count       native_count   
##  Min.   :    0   Min.   :    0   Min.   :    0.0   Min.   :  0.00  
##  1st Qu.: 1850   1st Qu.:  514   1st Qu.:  429.8   1st Qu.:  0.00  
##  Median : 5238   Median : 5168   Median : 1558.5   Median : 21.00  
##  Mean   : 9187   Mean   :10059   Mean   : 6848.5   Mean   : 54.08  
##  3rd Qu.:13503   3rd Qu.:15476   3rd Qu.: 8613.8   3rd Qu.: 70.75  
##  Max.   :60712   Max.   :69259   Max.   :69697.0   Max.   :601.00  
##                                                                    
##   asian_count      naturalized_citizen_count noncitizen_count uninsured_count
##  Min.   :    0.0   Min.   :    0             Min.   :    0    Min.   :    0  
##  1st Qu.:  279.5   1st Qu.: 2549             1st Qu.: 1450    1st Qu.:  643  
##  Median : 2151.0   Median : 6092             Median : 4157    Median : 1738  
##  Mean   : 4319.8   Mean   : 6625             Mean   : 4944    Mean   : 2325  
##  3rd Qu.: 5694.8   3rd Qu.: 9522             3rd Qu.: 7195    3rd Qu.: 3307  
##  Max.   :41314.0   Max.   :31847             Max.   :25008    Max.   :12182  
##                                                                              
##   school_count    eviction_count     store_count       sub_count     
##  Min.   : 1.000   Min.   :   1.00   Min.   :  1.00   Min.   : 1.000  
##  1st Qu.: 2.000   1st Qu.:  38.75   1st Qu.: 15.50   1st Qu.: 1.000  
##  Median : 5.000   Median : 141.50   Median : 40.00   Median : 1.000  
##  Mean   : 6.653   Mean   : 224.54   Mean   : 49.85   Mean   : 2.289  
##  3rd Qu.: 9.000   3rd Qu.: 358.00   3rd Qu.: 73.75   3rd Qu.: 3.000  
##  Max.   :31.000   Max.   :1126.00   Max.   :202.00   Max.   :17.000  
##                                                                      
##    bus_count      mean_ridership   perc_covered_by_transit          geometry  
##  Min.   :  1.00   Min.   :   273   Min.   :  0.00          MULTIPOLYGON :242  
##  1st Qu.: 30.00   1st Qu.:  4984   1st Qu.:  0.00          epsg:4269    :  0  
##  Median : 49.50   Median :  8006   Median : 26.40          +proj=long...:  0  
##  Mean   : 53.95   Mean   : 12051   Mean   : 41.86                             
##  3rd Qu.: 70.00   3rd Qu.: 14205   3rd Qu.: 64.17                             
##  Max.   :243.00   Max.   :109922   Max.   :254.34                             
##                   NA's   :119

Data Visuals

library(egg)
ggarrange(subway_loc, bus_loc, stops, bus_stops, ridership, access, ncol=2)

ggarrange(red, orange, yellow, green, teal, blue, purple, pink, brown, ncol=3)

ggarrange(white, black, latinx, asian, ncol=2)

Model Building

Models 1 & 2: No Grouping

We fit ungrouped poisson and negative binomial models below:

poisson_non_hierarchical <- stan_glm(
  white_count ~  total_pop + mean_income + mean_rent + 
    unemployment_count + sub_count + 
    perc_covered_by_transit + school_count + 
    store_count + bus_count + 
    eviction_count + uninsured_count,
  data = nyc_clean,
  family = poisson,
  chains = 4, iter = 500*2, seed = 84735, refresh = 0
)

negbin_non_hierarchical <- stan_glm(
  white_count ~  total_pop + mean_income + mean_rent + 
    unemployment_count + sub_count + 
    perc_covered_by_transit + school_count + 
    store_count + bus_count + 
    eviction_count + uninsured_count,
  data = nyc_clean,
  family = neg_binomial_2,
  chains = 4, iter = 500*2, seed = 84735, refresh = 0
)
check_1 <- pp_check(poisson_non_hierarchical) + 
  xlab("White Resident Count") +
  labs(title = "Poisson")+
  xlim(0,100000)+
  theme(plot.title =  element_text(face="bold", size=25, hjust=.5)) 
check_2 <- pp_check(negbin_non_hierarchical) + 
  xlab("White Resident Count") +
  labs(title = "Negative Binomial")+
  xlim(0,100000)+
  theme(plot.title =  element_text(face="bold", size=25, hjust=.5)) 

ggpubr::ggarrange(check_1, check_2, ncol=2, common.legend=TRUE, legend = "right")

nyc_predict_clean <- nyc_clean %>%
  na.omit()
set.seed(84735)

predictions_poisson <-  posterior_predict(
  poisson_non_hierarchical, newdata = nyc_predict_clean)

predictions_negbin <-  posterior_predict(
  poisson_non_hierarchical, newdata = nyc_predict_clean)

library(tidybayes)
library(bayesplot)

ppc_intervals(nyc_predict_clean$white_count, yrep = predictions_poisson,
              prob_outer = 0.95) +
  ggplot2::scale_x_continuous(
    labels = nyc_predict_clean$nta_id,
    breaks = 1:nrow(nyc_predict_clean)) +
  xaxis_text(angle = 90, hjust = 1) +
  theme_linedraw()+
  theme(panel.grid.major = element_line("transparent"),
        axis.title.y = element_blank(),
        axis.text.y = element_blank())+
  coord_flip()

ppc_intervals(nyc_predict_clean$white_count, yrep = predictions_negbin,
              prob_outer = 0.95) +
  ggplot2::scale_x_continuous(
    labels = nyc_predict_clean$nta_id,
    breaks = 1:nrow(nyc_predict_clean)) +
    xaxis_text(angle = 90,  hjust = 1) +
  theme_linedraw()+
  theme(panel.grid.major = element_line("transparent"),
        axis.title.y = element_blank(),
        axis.text.y = element_blank())+
  coord_flip()

Models 2 & 3: Hierarchy by Neighborhood

poisson_hierarchical <- stan_glmer(
 white_count ~  total_pop + mean_income + mean_rent + 
    unemployment_count + sub_count + 
    perc_covered_by_transit + school_count + 
    store_count + bus_count + 
    eviction_count + uninsured_count + (1 | nta_id),
  data = nyc_clean,
  family = poisson,
  chains = 4, iter = 500*2, seed = 84735, refresh = 0
)

negbin_hierarchical <- stan_glmer(
 white_count ~  total_pop + mean_income + mean_rent + 
    unemployment_count + sub_count + 
    perc_covered_by_transit + school_count + 
    store_count + bus_count + 
    eviction_count + uninsured_count + (1 | nta_id),
  data = nyc_clean,
  family = neg_binomial_2,
  chains = 4, iter = 500*2, seed = 84735, refresh = 0
)
check_1 <- pp_check(poisson_hierarchical) + 
  xlab("White Resident Count") +
  labs(title = "Poisson")+
  xlim(0,100000)+
  theme(plot.title =  element_text(face="bold", size=25, hjust=.5)) 
check_2 <- pp_check(negbin_hierarchical) + 
  xlab("White Resident Count") +
  xlim(0,100000)+
  labs(title = "Negative Binomial")+
  theme(plot.title =  element_text(face="bold", size=25, hjust=.5)) 

ggpubr::ggarrange(check_1, check_2, ncol=2, common.legend=TRUE, legend = "right")

library(kableExtra)
tidy(poisson_hierarchical, effects = "fixed", conf.int = TRUE, conf.level = 0.95) %>% 
  
  mutate(estimate= ifelse(term == "(Intercept)", exp(estimate), (exp(estimate)-1)*100), 
         conf.low= ifelse(term == "(Intercept)", exp(conf.low), (exp(conf.low)-1)*100), 
         conf.high = ifelse(term == "(Intercept)", exp(conf.high), (exp(conf.high)-1)*100))%>%
  filter(conf.low   > 0 & conf.high > 0 | conf.low  < 0 & conf.high < 0) %>%
  kable(align = "c", caption = "Hierarchicahl Poisson - Model Summary") %>% 
  kable_styling()
Hierarchicahl Poisson - Model Summary
term estimate std.error conf.low conf.high
(Intercept) 1717.1091298 0.3536747 885.2237713 3675.8332898
total_pop 0.0067707 0.0000086 0.0050549 0.0085170
mean_income 0.0010414 0.0000054 0.0000112 0.0021181
eviction_count -0.3087356 0.0005409 -0.4098675 -0.1987937
uninsured_count -0.0162629 0.0000601 -0.0274034 -0.0045796
tidy(negbin_hierarchical, effects = "fixed", conf.int = TRUE, conf.level = 0.95)%>% 
  
  mutate(estimate= ifelse(term == "(Intercept)", exp(estimate), (exp(estimate)-1)*100), 
         conf.low= ifelse(term == "(Intercept)", exp(conf.low), (exp(conf.low)-1)*100), 
         conf.high = ifelse(term == "(Intercept)", exp(conf.high), (exp(conf.high)-1)*100))%>%
  filter(conf.low   > 0 & conf.high > 0 | conf.low  < 0 & conf.high < 0) %>%
  kable(align = "c", caption = "Hierarchichal Negative Binomial - Model Summary") %>% 
  kable_styling()
Hierarchichal Negative Binomial - Model Summary
term estimate std.error conf.low conf.high
(Intercept) 2862.8732672 0.3439035 1459.0548198 5555.1716654
total_pop 0.0060187 0.0000086 0.0044313 0.0077559
mean_income 0.0010984 0.0000048 0.0001373 0.0021208
eviction_count -0.2514231 0.0004662 -0.3436878 -0.1586808
uninsured_count -0.0177119 0.0000511 -0.0279320 -0.0078445
nyc_predict_clean <- nyc_clean %>%
  na.omit()
set.seed(84735)

predictions_poisson <-  posterior_predict(
  poisson_hierarchical, newdata = nyc_predict_clean)

predictions_negbin <-  posterior_predict(
  negbin_hierarchical, newdata = nyc_predict_clean)

library(tidybayes)
library(bayesplot)

ppc_intervals(nyc_predict_clean$white_count, yrep = predictions_poisson,
              prob_outer = 0.95) +
  ggplot2::scale_x_continuous(
    labels = nyc_predict_clean$nta_id,
    breaks = 1:nrow(nyc_predict_clean)) +
  xaxis_text(angle = 90, hjust = 1) +
  theme_linedraw()+
  theme(panel.grid.major = element_line("transparent"),
        axis.title.y = element_blank(),
        axis.text.y = element_blank())+
  coord_flip() 

ppc_intervals(nyc_predict_clean$white_count, yrep = predictions_negbin,
              prob_outer = 0.95) +
  ggplot2::scale_x_continuous(
    labels = nyc_predict_clean$nta_id,
    breaks = 1:nrow(nyc_predict_clean)) +
    xaxis_text(angle = 90,  hjust = 1) +
  theme_linedraw()+
  theme(panel.grid.major = element_line("transparent"),
        axis.title.y = element_blank(),
        axis.text.y = element_blank())+
  coord_flip()